Exact Results for Propagators in the Geometrical Adhesion Model 
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The Geometrical Adhesion Model (GAM) that we described in previous papers provides a fully 
solved model for the nonlinear evolution of fields that mimics the cosmological evolution of pressure- 
less fluids. In this context we explore the expected late time properties of the cosmic propagators 
once halos have formed, in a regime beyond the domain of application of perturbation theories. 
Whereas propagators in Eulerian coordinates are closely related to the velocity field we show here 
that propagators defined in Lagrangian coordinates are intimately related to the halo mass func- 
tion. Exact results can be obtained in the ID case. In higher dimensions, the computations are 
more intricate because of the dependence of the propagators on the detailed shape of the underlying 
Lagrangian-space tessellations, that is, on the geometry of the regions that eventually collapse to 
form halos. We illustrate these results for both the ID and the 2D dynamics. In particular we 
give here the expected asymptotic behaviors obtained for power-law initial power spectra. These 
analytical results are compared with the results obtained with dedicated numerical simulations. 

PACS numbers: 98.80.-k, 98.80.Bp, 98.65.-r, 47.27.Gs 



I. INTRODUCTION 

Precision measurements and precision calculations of 
the statistical properties of the large-scale structure of 
the universe is becoming a key topic in observational 
and theoretical cosmology, in particular in the context 
of dark energy searches. Within the concordant model 
of cosmology (see e.g. [l[ and now strongly supported 
by the observations, see for instance Q), the dynamical 
growth of structure is a priori governed since z ~ 3000 
by collision-less dark matter components (although dark 
energy dominates the expansion at about z < 1), and 
structures are thought to emerge out of primordial nearly 
scale-invariant Gaussian metric perturbations. This gives 
rise to a hierarchical evolution, as increasingly larger 
scales turn nonlinear (very large scales being in the linear 
Gaussian regime and small scales in the highly nonlinear 
regime). Today, scales beyond ~ 10 Mpc are well de- 
scribed by linear theory while scales below ~ 1 Mpc are 
within the highly nonlinear regime. Providing theoretical 
insights into such a field is therefore a challenging idea. 

The highly nonlinear regime has proved very difficult 
to handle by analytical tools so far (see for instance [|[ ) , 
and one must resort to numerical simulations and phc- 
nomenological models (such as the halo model |4|) that 
involve some free parameters that are fitted to numerical 
results. 

The quasilinear regime on the other hand provides us 
with a priori robust and controlled predictions. In this 
regime indeed perturbation theory techniques can be em- 
ployed and many have already been developed. They are 
various extensions of the standard Perturbation Theory 
[3, such as the RPT 0-0], large- N expansions 0,1], the 
closure theory [l(J GjJ j the time renormalization equation 
[HI, etc., usually developed in Eulerian space but some- 
times in Lagrangian space as in [l3[ ■ The validity regime 



of such perturbation theories is however limited in na- 
ture, and in a non-controllable way, by the effects of shell 
crossings. Indeed, as small scale nonlinear regions form, 
fluid flows originating from different parts of the universe 
cross each-other leading to local multi-valued velocities, 
hence to vorticity and effective anisotropic pressure (see 
review Q). So far those effects have been observed to 
have limited effects on large scales, either through simple 
analytic investigation, as in flil . [l5| , or a posteriori from 
the successes of perturbation theory results. It remains 
that being able to explore the analytical properties of 
tools of interest in regimes were shell crossing is present 
is of crucial importance. 

The Geometrical Adhesion Model (GAM) provides for 
such an appealing framework [l|| [l7|. It is based on 
the Burgers equation (l8j in the inviscid limit which, be- 
fore shell crossing, reproduces the well-known Zel'dovich 
approximation [19J. At shell crossing however particles 
are prevented from crossing with the introduction of an 
infinitesimally small viscosity. The GAM is based on a 
subsequent prescription concerning the way matter flows 
within these critical regions. This model is examined in 
[20l l2l| where it is argued that it provides an attractive 
toy model for the cosmic matter distribution. In this pa- 
per we aim at taking advantage of this fact to explore 
statistical indicators in their full complexity. 

The quantities we are more particularly interested in 
are the so-called propagators. Indeed, while standard 
perturbation theory only involves many-body density 
and velocity correlation functions, or polyspectra, most 
resummation schemes also involve "propagators" or "re- 
sponse functions" . These unequal-time quantities de- 
scribe the evolution with time of small fluctuations, and 
often appear at intermediate stages in these resumma- 
tion procedures. Such quantities appear for instance in 
the "large- N v expansions developed in d,H,[ll|, the "clo- 
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sure approximation" of [l(| [H| and also in the "RPT" 
approach of [5H2( where they were initially introduced. In 
addition to encapsulating diagrammatic resummations, 
they also have a well-defined physical meaning. Prop- 
agators are defined, as is usual in statistical physics, 
as two-time response functions of the system. For in- 
stance the density propagator in real space coordinates, 
G<k(x2, ^2! Xi, ii), is the ensemble average, over all inter- 
vening modes, of the functional derivative 



t 2 >h: G5c(x 2 ,t 2 ;xi,ti) 



T>6(x 2 ,t 2 ) 
2>C(xi,ii) 



(1) 



where £ is an infinitesimal external perturbation to the 
density contrast, or velocity field, applied at time t\. This 
gives the "linear response" of the system, with respect 
to the density contrast at position x 2 and time t 2 , to 
an external perturbations applied at position xi at the 
earlier time t\ [49]. In this article, we shall focus on the 
propagator from initial-time fluctuations, t\ — ¥ 0. Then, 
this is the response to the initial conditions themselves, 
which can be defined through the initial velocity potential 
■00 or the linear density contrast Slo- Thus, as in [23j j . 
we define the density propagator, 



t>0: G 4 (x,t;q) = 



^Lo(q) 



(2) 



It gives the change in the density contrast at position 
x and time t induced by an infinitesimal change to the 
initial (or linear) density contrast at position q. 

One of the aims of this paper is to derive the prop- 
erties of propagators in as much details as possible and 
compare the results with numerical experiments. In par- 
ticular we want to see how these functions behave, and 
what they are sensitive to, when a full nonlinear evo- 
lution of the fields is taken into account. There are a 
priori two means of describing the large-scale structure 
dynamics and properties; one is based on the use of the 
Eulcrian coordinates, this is a natural choice since it is 
directly related to observations, the other is based on the 
Lagrangian coordinates. As we shall see below, and as 
pointed out in [22J, |24j , the Eulerian response functions 
and propagators are dominated by a "sweeping effect", 
that is, the collective transport of small-scale density fluc- 
tuations by the long wavelengths of the velocity field . 
This holds for both the gravitational and the Burgers 
dynamics, and it means that the properties of these two- 
time functions are not a very good probe of the prop- 
erties of the large-scale structures, observed at a given 
time. This shortcoming led to the study in [2^, [2J] of 
Lagrangian-space propagators. Indeed, in a Lagrangian 
framework, where one follows the motion of particles, the 
impact of uniform translations is automatically removed 
(because they do not change the system as viewed from a 
particle). Then, propagators or correlation functions au- 
tomatically go beyond the "sweeping effect" and directly 
probe the deformation of the density field (e.g., tidal ef- 
fects). This suggests that Lagrangian propagators should 



generally provide more sensitive probes of the matter dis- 
tribution. This is one of the motivations for the study 
presented here. In the context of the GAM it is indeed 
possible to relate propagators in Lagrangian coordinates 
to geometrical properties of the late-time field. 

The plan of the paper is the following. In Sect. [TT| we 
present the equations governing the Burgers dynamics 
and the geometrical adhesion model. Analytical results 
regarding this model are presented in the following sec- 
tion IIII1 In Sect. IIVI we present numerical results that 
illustrate our anaytical findings. 



II. BURGERS DYNAMICS AND 
GEOMETRICAL ADHESION MODEL 

A. Equations of motion and geometrical 
constructions 

We first recall in this section the definition of the 
"geometrical adhesion model", GAM, described in more 
details in a previous paper [20|. This model coincides 
with the well-known Zel'dovich dynamics [l9| before shell 
crossing and beyond shell crossing it is built from 

the d-dimensional Burgers equation [l8[ for the velocity 
field u(x, t), 



dtu + (u • V)u = i/Au, 



(3) 



where we consider the inviscid limit. As is well known, 
for curlfree initial velocity fields the nonlinear Burgers 
equation fiSb c an be solved through the Hopf-Cole trans- 
formation [2^, [26[, by making the change of variable 
ip(x,t) = 2v lnH(x, t), where ijj(x,t) is the velocity po- 
tential defined by 



u(x,t) = -VV>. 



(4) 



This yields the linear heat equation for 3(x, t), which 
leads to the solution 



^(x,t) = 2v In 



dq 



exp 



V'o(q) 

2v 



|x-q|- 



Then, in the inviscid limit v — > + , a steepest-descent 
method gives pj, H3 



^(x,t) 



sup 
q 



(q) - 



2t 



(6) 



If there is no shock, the maximum in ^ is reached at a 
unique point q(x, t), which is the Lagrangian coordinate 
of the particle that is located at the Eulerian position 
x at time t (hereafter, we note by the letter q the La- 
grangian coordinates, i.e. the initial positions at t = of 
particles, and by the letter x the Eulerian coordinates at 
any time t > 0) . If there are several degenerate solutions 
to ©, we have a shock at position x and the velocity is 
discontinuous. 
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Thus, a key property of the nonlinear equation ([3]) 
is that we know its explicit solution © at any time t. 
Therefore, we can build the velocity field at any time 
t from the maximization ([S]), which also corresponds to 
a Legendre transform, without solving the dynamics at 
intermediate times. In a cosmological context we are 
more particularly interested in the matter distribution; 
it is therefore very useful to extend this property to the 
density field. As explained in detail in [20] , this is possi- 
ble provided we use a specific continuity equation, which 
differs from the usual continuity equation by a peculiar 
diffusive term, proportional to v, that vanishes outside of 
shocks in the inviscid limit. More precisely, one obtains 
in this case the matter distribution from the Lagrangian 
map (i.e. the displacement field), qn>x, which is defined 
as follows. 

Defining the "linear" Lagrangian potential ip l (q, t) by 



Then, the densityfield is determined by the conserva- 
tion of matter [H IM HI, 



p(x,t)dx = p dq, 



which reads as 
P(x) 



Po 



dot 



9q 

Ox 



dot 



dx 
<9q 



(14) 



(15) 



Here we used the fact that both determinants are posi- 
tive, thanks to the convexity of i?(x) and y(q), and we 
assumed a uniform initial density p^ at t = 0. Thus, 
the "Lagrangian-Eulerian" mapping q f-> x of Eq. lfTTj) 
and the density field (fT5j) define the "geometrical adhe- 
sion model" that we study in this article, where both the 
velocity and density fields can be obtained at any time 
from the initial fields by geometrical constructions. 



^(q,*) = ^--*Mq) 



(7) 



so that in the linear regime the Lagrangian map is given 
by 



Xi(q,t) = -q^- = q + <u (q), 



and introducing the function 



H(x,t) 



ttp{x,t), 



(8) 



(9) 



one can see that the maximumj|6|) can be written as the 
Legendre transform (see note |5l|) 



H(x,t) 



sup 
q 



xq 



lor 

2 



#o(q) 



(10) 

Therefore, Eq. (|10p yields the inverse Lagrangian map, 
x n- q, q(x, t) being the point where the maximum in 
Eq.© or (|10[) is reached. This is only a rewritting of 
the Hopf-Cole solution ©, but this is not sufficient to 
fully define the density field in dimension greater than 
one. Indeed, because of shocks there is no unique way 
to invert the mapping x H> q. Then, the "geometrical 
adhesion model" [13, [2(| consists in choosing the direct 
Lagrangian mapping, q i— > x, as the Legendre conjugate 
of x i Y q, 



, , dH , , dip 
q(x,i ) = -5- , x(q,i) = — , 
Ox Oq 

where the potential (p is given by 



(11) 



<p(q, t) = £ q [iJ(x, *)] = sup [q • x H{x, t)} . (12) 

X 

From standard properties of the Legendre transform, this 
only gives back the linear Lagrangian potential <p l (q, t) 
of Eqs.© and (TTU)) if the latter is convex, and in the 
general case it gives its convex hull, 



p> = conv(<pi). 



(13) 



B. Initial conditions and self-similarity 

As in [2~J | , we consider Gaussian random initial condi- 
tions, which are isotropic and homogeneous, with power- 
law initial power spectra. They can be defined in terms 
of the initial velocity potential, ^o(x), or equivalently in 
terms of the linear density contrast, <5l(x, t), which is 
the usual approach in cosmology. Thus, introducing the 
density contrast, 



5(x,t) 



p{x,t) - po 



Po 



(16) 



one can see that in the linear regime [52| , where the parti- 
cles follow the trajectories (j8|) that still coincide with the 
Zel'dovich dynamics [l^] , the density contrast behaves at 
linear order as 

5 L (x,t) =tS L0 (x), with (5 L0 = -V-u = AV>o- (17) 
Going to Fourier space, with the normalization, 



S L (x,t) 



dk e ik x <J L (M), 



(18) 



the linear density field 5l is taken as Gaussian, homo- 
geneous, and isotropic. Then, it is fully described by its 
power spectrum Pg L (k,t), 

(Sl)=0, C5 L (k 1 )~6 L (k 2 ))=5 D (k 1 +k 2 )P SL (k 1 ), (19) 

which we choose of the power-law form 

D 



3<n<l: Pg L (k,t) = 



{2ir) d 



^2 



(20) 



where D and n are the amplitude and slope parameters. 

In the inviscid limit, one can check [U |3l| that 
the power-law initial conditions (f2U|) give rise to a self- 
similar dynamics [53| : a rescaling of time is statistically 
equivalent to a rescaling of distances, as 



A>0: i->At, x -> A 2/(n+3) x. 



(21) 



4 



Thus, the system displays a hierarchical evolution and 
the only characteristic scale at a given time t is the so- 
called integral scale of turbulence, L(t), which is gener- 
ated by the Burgers dynamics and grows with time as in 
(|2~T1). Hereafter we choose the normalization 



L(t) = (2Dr) 



2U/(™+3) 



(22) 



where the constant D was defined in Eq.(f2"U|). This scale 
measures the typical distance between shocks, and it 
separates the large-scale quasi-linear regime, where the 
density power spectrum keeps its initial power-law form, 
(1201) , from the small-scale nonlinear regime, which is gov- 
erned by shocks and pointlike masses, where the density 
power spectrum reaches the universal white-noise behav- 
ior (i.e. Ps(k,t) has a finite limit for k ^> 1/L(t)). De- 
tailed numerical studies of the cluster mass functions and 
density fields generated in ID and 2D are presented in 

[ilEj. 

III. ANALYTICAL RESULTS 
A. Definitions 



Multiplying by the operator Cq we recover Eq. ((23)) . This 
result obtained here for Gaussian initial conditions can 
obviously be extended to non-Gaussian initial conditions 
(but then the integration by parts gives rise to a more 
complicated functional form of the cross-correlation with 
the initial fields). What this simple calculation shows 
is that the identities (|23|) and (|27[) are quite general. 
They apply to non-analytic functionals J- X _t and should 
be valid as soon as propagators actually exist. In particu- 
lar, these relations remain valid after shell crossing, where 
perturbative expansions break down. This is confirmed 
in App. [C] where we do not detect any deviation from 
the identities (|23|) and (f27|) far in the nonperturbative 
regime. The full validity of this identity was overlooked 
in early papers but was successfully checked in numeri- 
cal simulations in Q. This is a useful property since it 
is usually more convenient to measure cross-correlations 
such as (<5(x, t)5Lo(<l)) m numerical simulations whereas, 
in the case we consider here, the expression @ is better 
suited to analytical investigations. 



2. Eulerian quantities 



1. Propagators and correlators, general properties 

As noticed in [f| from the perturbative expansion of the 
nonlinear density contrast <5(x, t) over powers of the lin- 
ear growing mode Slo, for Gaussian initial conditions the 
response function ([2]) is related to the cross-correlation 
&(x,i;q) by 



(J(x,i)fco(q)) 

dq'G 6 (x,wx<Mq')Mq)>-(23) 



&(x,t;q) 



Following the approach used in [22j , this relation can ac- 
tually be obtained by a simple integration by parts and 
does not explicitly rely on perturbative expansions. In- 
deed, writing the nonlinear density contrast as a func- 
tional of the initial condition Slq, 



S(x,t) = J" x ,t[(5Lo(q)] 
the response function ([2]) writes as 



G 4 (x,t;q) = / V5 



5lo e 



(24) 



, (25) 



where we did not write an irrelevant normalization con- 
stant and C (qi,q 2 ) = (<5z,o(qi)<Wq2)) is the two- 
point correlation of the initial Gaussian field. Integrating 
Eg. (f2"5|) by parts gives 



G s (x,i;q) = Jv8 L0 e-^ SL °- c ° , t [6i 
J dq'C - 1 (q,q')-^o(q') 



(26) 



= J dq'C - 1 (q,q')-(<5(x,t) ( 5 L o(q')).(27) 



In the following we will be interested in quantities de- 
fined in both Eulerian coordinates and Lagrangian co- 
ordinates. In practice, since the system is statistically 
homogeneous and isotropic, it is convenient to work in 
Fourier space, and we define the Fourier-space counter- 
part of the response Gs of Eq.Q by 

^ 0: U:f k ;,t )=S D QL-V)G 6 (k,t), (28) 
\X>dio(k'J / 

which is related to the real-space response by 



G s {k,t) = / dxe- ikx G s (x,t), 



(29) 



where we used Gg(x, i;q) = Gg(\x — q|,t). The relation 
reads in Fourier space as 



S(k,t)S L0 (k')) =5 D (k + k')G s (k,t)P L0 (k). (30) 



Here we focused on the density propagator or response 
function, but we can also consider the velocity potential 
ip, such as 



G^(x,i;q) 



£ty>o(q 

and mixed statistics, such as 

_ / x>a(x,t) \ 

G ^-\PMqT/' 
and their respective Fourier-space counterparts. 



(31) 



(32) 
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3. Lagrangian quantities 

Similar quantities can be defined in a Lagrangian 
framework. Following [24j we first introduce the diver- 
gence k of the displacement field, 



n(q,t)=d- — 



-V q -[x(cbi)-q], (33) 



where x(q, t) is the position at time t of the particle of 
Lagrangian coordinate q. From the conservation of mat- 
ter, Eq. (fT5|) , we can see that at linear order k is also the 
density contrast, 

K L (q, t) = S L (q, t) = t Klo (q) with k lo (q) = 5 L0 (q) . 

(34) 

It is clear from the definition (|33|> that a uniform trans- 
lation does not change the value of k, so that it is not 
sensitive to the "sweeping effect" discussed in the intro- 
duction [22, 24]. Then, as above, we define the associated 
Lagrangian propagators, 



t>0: G K ^(q,t;q ) 



P«(q,f) 
X'V'Lo(qo) 



and 



i>0: G K (q,i;q ) 



Vn(q,t) 
Vn L0 (q a ) 

This again reads in Fourier space as 



(35) 



(36) 



/ T>k(k,t) 
t > : < ■ V ' ; 



and 



* > 



\X>V>Lo(k') 

/ D«(k,i) 
\X>KLo(k') 



= Mk-k')G KV ,(M), (37) 



= 5z?(k-k')G K (fc,t), (38) 



Using the property ^lo(^) = ^Lo(fc) = —ki/)o(k), from 
Eqs. (|34II7p . we also have the relation 

G K (M) = -fc- 2 G K ^(fc,t). (39) 

Similarly to the Eulerian case we also have the relation 

(k(k, i)KLo(k')) = S D (k + k') G K (fc, t) P L0 (k). (40) 

Thanks to the explicit Hopf-Cole solution ([6]) and its ge- 
ometrical interpretation in terms of first-contact parabo- 
las [lH H3] , very fruitful insights on the behavior of the 
propagators can be obtained. Exact results for the one- 
dimensional case have been recently derived in [23j , which 
we will briefly recall in the course of the general calcula- 
tion. 

B. Eulerian response functions 

To take advantage of the Hopf-Cole solution ^ it is 
convenient to first consider the response function G^, as 
defined in Eq.flHl)- 



In any dimension we can use the Hopf-Cole solution 
dHJ) to obtain the propagator G^. Thus, as in the ID 
case [23j |. taking the functional derivative of Eq.© and 
next taking the inviscid limit we obtain 



G v ,(x,t;q ) = (<Mq(x,t) -qo]) =px(q ,t) 
= ^P(u,t), 



with 



qo 



t 



(41) 



(42) 



Here p x (qo,i) is the probability distribution function of 
the Lagrangian coordinate qo of the particle that is lo- 
cated at position x at time t, and p(u, t) is the probability 
distribution function of the velocity u at position x and 
time t (it does not depend on x because of statistical 
homogeneity) . 

Therefore, the Eulerian propagator is always given 
by the one-point velocity probability distribution, irre- 
spectively of the space dimension. As such, it is still 
governed by the "sweeping effect" mentioned in the in- 
troduction, and for convergent initial power spectra it 
obeys in the weakly nonlinear regime behavior, 



~(RPT 



(M) 



cr 1 (*>*;<i) 



, -t 2 k 2 a 2 12 

1 



{V2^ta UQ ) d 



-q| a /(« a O j (43) 



where a^ o = (|uo| 2 )/<i is the variance of the initial veloc- 
ity along any given direction. The behavior (|43[) , which 
is more general than the adhesion model, also applies to 
the gravitational dynamics in general. Indeed, it only 
relies on the transport of particles by long-wavelength 
modes of the velocity field (see 0, |32| for an account of 
this effect in a perturbative approach), in other words, on 
the Galilean invariance further discussed in Sect. MI C 21 
below. It is also at the basis of the RPT resummation 
scheme introduced in cosmology 

As discussed in [23| for d = 1, and remains valid in 
higher dimensions, for the power-law initial conditions 
(|20p the weakly nonlinear regime (|4"3")l does not exist. 

For — 3 < n < — 1, the variance a\ diverges because 
of low-A; modes. This means that as we increase the size 
of the system (since in practice, for instance in numeri- 
cal simulations, we consider finite-size systems and even- 
tually take the infinite-size limit) the contribution from 
long wavelengths keeps increasing and particles are trans- 
ported over increasingly large distances. This yields a 
divergent "sweeping effect" and the propagators vanish 
as soon as t > 0. However, this is not a genuine loss of 
memory since this divergence is due to almost uniform 
random translations and thanks to Galilean invariance 
the structures of the density field are not affected. In 
particular, as we recall in Sect. iHXCl below. Lagrangian 
propagators remain finite. 

For — 1 < n < 1, the variance cr 2 diverges because of 
high-fc modes. However, as soon as t > this ultraviolet 
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divergence is regularized by the dynamics, more precisely 
by the nonperturbative "sticking" of the particles at colli- 
sion. This means that the one-point velocity probability 
distribution p(u,t), whence the response functions, are 
finite and well-defined. However, they are governed by 
nonperturbative effects and do not obey the Gaussian be- 
havior (|4"3")l . For d — 1 and in the case n = one can actu- 
ally obtain its exact expression [23|, [3^, 34j , and this gives 
in particular the exponential decays G ~ e~( x ~ q ' I 1 ' and 
G ~ e~ tk 1 at large \x — q\ and large k. For generic n 
we only know the large-separation tail [3l|, HH, [36| 



1 < n < 1, \x— q\ — » oo : G^(x,t;q) ~ e" 



.( 44 i 

which depends on n, contrary to the weakly nonlinear 
regime behavior (|4"3")l . For d > 1 we no longer have exact 
expressions for p(u,t), even for rt = 0, but in the large- 
separation limit we expect from rare-event analysis [31] 
the tail 



-1 < n < 1, |x— q| — > oo : G^,(x,t;q) ~ e" 



|x~q|<"+ 3 >/« 2 

(45)' 
as in flU}. 

The probability distributions p(qo, i) and p(u, i) being 
well defined and normalized to unity for — 1 < n < 1, we 
obtain from Eq. (|4"Tj) the low-fc limit 



Gv(M) = 1. 



(46) 



This agrees with the linear regime prediction but this 
exact result is more general since it follows from the non- 
perturbative identity (|41l) . Thus, it remains valid for the 
scale-invariant initial conditions with —1 < n < 1, where 
there is no true weakly nonlinear regime in the sense of 
(|4"3")l because linear velocity fields are divergent. How- 
ever, particles are not transported arbitrarily far away 
but on distances that scale with the characteristic length 
L(t) defined in (f2"2")) . Thus, even though the "confining" 
process is strongly nonperturbative, due to the collisions 
and sticking of particles, initial density fluctuations seen 
on scales much larger than L(t) are not erased and evolve 
as in linear theory (in a weak sense, that is, if we con- 
sider smooth windows such as a Gaussian rather than a 
top-hat). In terms of Burgers turbulence this is related 
to the principle of "permanence of large eddies" [3(| ■ 

For d = 1 the density contrast 5{x, t) associated with 
the matter distribution (IT51) also reads as 



S(x,t)=t 



d 2 *p 
dx 2 



(47) 



Then, from Eq. (|4"Tj) one obtains (e.g., by going to Fourier 
space) 



Gg(x, t; q) = p(u, t), with u 



x — q 
t 



(48) 



In agreement with physical arguments and perturbative 
analysis [H, [24[ , the results flU]) and (gSJ) explicitly show 
that Eulerian response functions are dominated by the 



"sweeping effect" , that is, by the transport of particles 
by the underlying velocity field. Since only the one-point 
velocity distribution appears in (|48|) . these Eulerian prop- 
agators are not very sensitive probes of the density field 
(which depends on relative motions). 

For d > 2 the density contrast is no longer related to 
the velocity potential -0 by a linear relationship. From 
Eas. plTg)) we obtain 



dxj 1,3 dxidxj 



and Eq.(fl5l) yields 



1 + 5(x,t) = dct 



(49) 



(50) 



which contains terms up to power d over ip. The Eule- 
rian density response function must still be sensitive to 
the "sweeping effect" and dominated by the properties of 
the velocity field, in agreement with the response func- 
tion of the velocity potential ip itself. However, because of 
the nonlinear dependence (|50p the response function G$ 
is no longer merely proportional to the one-point veloc- 
ity probability distribution and the relationship is more 
complex, and depends on the dimension d. For power-law 
initial conditions it should again vanish as soon as t > 
for —3 < n < — 1, and be finite for — 1 < n < 1, with the 
asymptotic behavior (|4"5|) at very large separations. 



C. Lagrangian response functions 

1. General formalism 

As seen in [2(1 [37| , for the power-law initial conditions 
that we consider here, all the matter clusters into a set of 
point masses that builds a Voronoi-like tessellation in Eu- 
lerian space and a dual triangulation in Lagrangian space. 
This Lagrangian-space tessellation is built from segments 
in ID, triangles in 2D, tetrahedra in 3D, and higher-order 
simplices in higher dimensions |28j . A sketch of such a 
partition in 2D is to be found on Fig. [TJ In this picture 
one halo is associated with each Lagrangian triangle and 
it is formed from all the particles that were initially in 
this corresponding triangle. As a result the mass distri- 
bution is nothing but the volume distribution of those 
triangular cells. Moreover, the position of the halos in 
Eulerian space depends on the potential heights Lp(qf^) 
on each of the summits of the triangle. As time 
evolves the relative heights of the summits change leading 
to halo motions and eventually to halo mergings. These 
processes were described in details in [201 ] . 

The final matter distribution can be characterized by 
the density field or cquivalcntly by the displacement field, 
a quantity that better suits the Lagrangian description. 
What this construction tells us is that within each La- 
grangian cell V^") the Lagrangian mapping x(q) = x( Q ) 
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FIG. 1: In a triangulation, the matter in each triangle moves 
toward a single halo. Infinitesimal variation of the potential 
at position q i (<») induces an infinitesimal change in the dis- 
placement field within the triangles q^a) is a summit of. This 
change can entirely be described by a change in the positions 
of the corresponding halos. 

is constant, as all particles that originate from this re- 
gion belong to a single point-like cluster at position x( Q ) 
at time t. And because the displacement is potential it is 
entirely determined by its divergence k. It is clear here 
that the divergence within each cell is constant. Indeed 
from the definition ([3"3"| this yields «(q, t) = d within 
each of these Lagrangian cells, with the addition of a 
Dirac term on the boundaries of these cells associated 
with the jump from to x' Q ' as one goes from cell 
V( a) to cell V (Q,) . This more formally reads as 

«(q) = d - f d d - x s S D [q - qjj(s)] |n ■ Ax| (51) 

where qs(s) is the (d — l)-dimensional manifold E built 
by the boundaries of the Lagrangian cells (i.e., all trian- 
gle sides in 2D), parameterized by a coordinate s, d d_1 s 
is the natural surface element on S embedded in d dimen- 
sions, n is the unit normal vector to this manifold and 
Ax = x^" ) — x( Q ) is the separation vector between the 
Eulerian positions x^ and x^ Q ) of the two mass clus- 
ters associated with the two Lagrangian cells (a) and (a') 
that are on either side of E. In ID the integral (f5"Tj) sim- 
plifies to the sum over all endpoints qi of the Lagrangian 
intervals (Zi+i], associated with the shocks of Eulerian 
position Xi, as n(q) = 1 — £\ S D (q - <fc) fa - Xi-i) [11]. 

It is convenient to decompose the integral (|5"Tj) over 
all cells V^ a \ The absolute value |n • Ax| arose from 
the convexity of the mapping x(q), that is, the fact that 
<p(q) is convex in Eq. fllip . This ensures that dxijdqi > 
along any direction i [TaHMIH^! so that the contribution 
(n • Ax) must be taken positive. This can also be written 
as 

-|n.Ax|=n^.xW+n^.x(«'), (52) 



where and are the unit normal vectors to the 

surface E' a ' a ) that point outward from the neighboring 
cells V< Q > and V {a ">. Then, Eq.® reads as 




ds-x^ 5 D [q-q s(a )(s)], (53) 



where we note ds = n out ds. 

Next, to estimate the Lagrangian response functions 
Q35I36I) we must consider the variation of Eq. (|55|) for in- 
finitesimal changes of the initial conditions. The mani- 
fold E is set by the convex hull construction (|13p . For 
instance, in ID the boundaries qi are the contact points 
between the linear Lagrangian potential Ph{q) and its 
convex envelope (p, whereas in 2D the triangles of the 
Lagrangian-space tessellation are the triangular facets of 
the convex hull p (and the triangle summits q^ are 
again the contact points between ipi and ip, for each facet 
(a) and i = 0,1,2 indexes the three summits). Then, in 
the nondegenerate case (e.g., in 2D a planar facet only 
makes contact with three points) , which has a unit prob- 
ability, an infinitesimal change of the initial conditions, 
whence of PL{q), does not modify the set {q^ } of the 

contact points, but only their heights p(q[ a ^). This is 
because for the power-law initial conditions (|20|) the po- 
tentials V'o an d Pl have no finite second-derivatives and 
the contact points are isolated infinitesimally-thin spikes, 
see [2(J. As a result, infinitesimal variations of the ini- 
tial conditions do not change the manifold E neither the 
Lagrangian cells V*-"-*; they only change the Eulerian po- 
sitions x(") of the mass clusters associated with each La- 
grangian cell. This is the core property which determines 
how Lagrangian response functions and halo mass func- 
tions are related. We then have 

X>K(q) f . £>x< Q ) 

^o(q') f-fJxi") ^Vo(q') 

(a) 

Following the definition (|35|) the Lagrangian propaga- 
tor is then given by the statistical average of Eg. ([5^]). 
Because the system is statistically homogeneous and 
isotropic, G rel /,(q, q') only depends on |Aq| with Aq = 
q q' , and we obtain 

G„^(Aq,t) = / dmn(m)( [ ds ■ — — — ) 
Jo \Js ^o(q(s)-Aq)/ m 

where n(m)dm is the mean number of Lagrangian cells 
of mass m (or of point-like clusters of mass m in Eulerian 
space), per unit volume. Here (..) m is the statistical av- 
erage over cells of fixed mass to. In ID this is immediate 
since Lagrangian intervals are fully defined by their mass 
(i.e. their length), but in higher dimensions we must av- 
erage over the distribution of shapes and angles at fixed 
mass to. We will show explicit examples of such averag- 
ing in the following. 

As recalled above, the Eulerian position x of the mass 
cluster only depends on the values ipo-i of the velocity 
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potential at the (d + 1) summits that define the La- 
grangian cell V, hence we can write the functional deriva- 
tive in Eq. (f55|) as 

T>x <9x . , 

^^ = g^ Mq - q '°- (56) 

By symmetry the contributions of all (d + 1) summits 
are identical, after we average over shapes and angles, 
and choosing the origin of coordinates on summit qo we 
obtain 

/>OC 

G K tP ( Aq, t) = (d+1) / dmn(m) 
Jo 

X (/ E ds '^ fo(< > (s) - A< > ) ) m ' (67) 

where qo = 0. This yields in Fourier space 

poo 

G^{k,t) = -tk 2 dmn(m)l d W(k;m), (58) 
Jo 

where I is the typical size of cells of mass m, defined as 
m = p t d , (59) 
and the dimensionless kernel W{k; m) reads as 

(60) 

where again (..) m is the statistical average over cells of 
fixed mass m. Then, from Eq.(f3l?f the propagator G K 
writes as 

/>oo 

G K (k,t)=t dmn(m)e d W{k;m). (61) 
J o 

2. Large scale behavior 

We first explore the behavior of W{k;m) for k — > 
and show that it correctly reproduces what one expects 
from the linear theory. For that we simply use Eq. ([60|) 
and Taylor expand it at k 0. Let us first note that the 
term in 1/k vanishes by symmetry, because the factors 
k ■ [q(s) — qi) average to zero as we integrate over the 
global orientation of the Lagrangian cell V. We then get 



W(k;m) = 




+0(k 2 ), (62) 

where we replaced the factor (d+ 1) by the sum over the 
derivatives with respect to the (d + 1) summits, going 



back to Ea. (l56|) . and where in the second term k = k/|k| 
is a unit vector along an arbitrary direction. It can then 
be observed that the sum in the first term vanishes since 
the position x of the mass cluster is not modified by 
a uniform shift of the velocity potential. We can also 
observe that, using Gauss' theorem, the integral on the 
surface E in the second term can be written in terms of 
an integral on the volume V. Using again the property 
J2i dx/di/jo-i = 0, this yields a volume factor |V| = l d so 
that 

^(fc;m) = ^/^(k- qi )(k-^-)\ +0(k 2 ). (63) 

This sum is in fact constrained by the Galilean invariancc 
of the dynamics. Indeed, adding a uniform initial velocity 
vo, Uo(q) — > uo(q) +vo, also corresponds to the changes 

x(q, t) x(q, t) + v t, Vo(q) -> ^o(q) - v • q (64) 

and does not modify the structure of the Lagrangian- 
space tessellation. This implies that an infinitesimal uni- 
form velocity perturbation vo leads to the shift of cluster 
position 

v i = Ax = V— — Aipo-i = - V ^- — (v -qi), (65) 
and taking the scalar product with Vo yields 

|vo| a = ^$>°-*)( v °-^)- (66) 

i 

This equation holds for Vo of any direction and any length 
(since both sides scale linearly with |vo| 2 ), and taking 
| v | = 1 we obtain that the average in (|63p writes as 
W(0;m) — (l) m , whence 

W(k;m) = 1 + C(fc 2 ), (67) 

independently of the statistical properties of the 
Lagrangian-space tessellation. 

Finally, since all the matter is contained in the mass 
clusters, the mass function obeys the normalization 

poo 

I dmn(m) — = 1, (68) 
Jo Po 

and using Eqs. (|59l61[) we obtain 

G K {0,t)=t, (69) 

which agrees with the linear regime prediction associated 
with Eg. (IMl) . i.e. the displacement field follows the linear 
theory at large enough scale. 

It is interesting to make the connection between the 
exact result (|69p , which does not rely on perturbative ar- 
guments, and the comparison between Eulerian and La- 
grangian propagators. The Galilean transformation (|64[) 
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merely states how particles are transported over a dis- 
tance v t by long-wavelength modes of the velocity field. 
In the Eulerian framework, this "sweeping effect" leads 
to the Gaussian decay (|4"3")l after we take the statistical 
average over these long-wavelength modes. In contrast, 
in the Lagrangian framework the analysis above explic- 
itly shows that this Galilean law only sets the normal- 
ization at k = of the propagator, and its values at 
k > probe how trajectories depend on the curvature 
and higher order derivatives of the velocity potential, as- 
sociated with higher orders of the expansion over k of 
the expression (|60l) . Moreover, in the Lagrangian frame- 
work the Galilean transformation (|64|) appears through 
relative displacements of well-separated regions instead 
of absolute displacements [HJ]. 



3. Scaling laws 

In arbitrary dimensions computing the averages (|5T[) 
or (|60p can lead to lengthy expressions. However, before 
we explicitly consider the ID and 2D cases, we give some 
general scaling arguments. We first need to evaluate the 
variation of the halo position x for infinitesimal variation 
of the potential at a summit position . By construction, 
in the (d+1) space {q, <p} the convex hull f(q) of Eq. (fT3")) 
is an hyperplane over the region V, defined by the (d+1) 
points {qi,v?L;i}, and its equation reads as 



dot 



qi qo 



qd - qo q qo 



= 0. 



Thus, over this cell <£>(q) takes the form 

¥>(q) = </?i;o + (v^l) • (q - qo), 



(70) 



(71) 



where we note (V</?l) a constant vector that is linear over 
the potentials {</>L;0, ■ •, fL-d} and that scales as the ratio 
(6<pl)/£, where £ and (5<pl) are the typical differences of 
position and height between the (d+1) summits (£ was 
defined in Eq.([59"])). From Eq. flTTj) the position x of the 
associated mass cluster is given by x = (V</2l), and from 
Eq.([7]) we obtain the scaling 



9x 



(72) 



Substituting into Eq. ([60|) we check that the kernel 
W(k; m) is dimensionless and behaves as 



(73) 



W(k-m) ~ W (k£), 



in terms of a scaling kernel Wq. We have seen in Eq. (|67j) 
that W(0;m) = 1, whence W (0) = 1. On the other 
hand, at high k we can expect the exponential factor in 
Eq. flrJUl) to cut large distances beyond s ~ 1/fc, which 
gives a power-law decay, 



At low masses, numerical simulations and heuristic argu- 
ments [13, HH H3] suggest that the cluster mass function 
obeys the power-law tail 



3 < n < 1, m < poL(t) d : n(m, t) ~t 



i)/2 j 

which has only been proved rigorously in ID for the 
Brownian case n = — 2 (4p| - |42l | and for the white- noise 
case n = (13, Hill- Then, the cutoff $J4$ gives rise to 
two different behaviors. For small values of n the integral 
(|6T|) is dominated by small masses m, with m oc i d and 
I — 1/k, because of the cutoff ((74]) of Wo(k£), and we 
obtain the scaling 

n<--l, fc>L(t)- 1 : G K (k,t)~t(kL)- d(n+3 ^ 2 , 
d 

(76) 

whereas for large values of n the cutoff (|T4"1) is too shallow 
and the integral is dominated by large masses (set by the 
cutoff of the mass function), and lengths I ~ L(t) of the 
order of the typical nonlinear scale ([22]) , which yields 

rc>~-l, k^>L(ty l : G K (k,t) ~t(kL)- d -\ (77) 



The scaling laws (|76|77j) have been derived assuming 
that small-mass Lagrangian cells are characterized by a 
single scale £, which leads to Eq. ([73|) . As we discuss in 
Sect. [TO C 51 below, in dimensions greater than one small 
cells may be characterized by several lenght scales, which 
behave as power laws over mass with different exponents. 
This would violate the scaling laws (|76I77[) and give rise 
to new exponents. 



4- ID case 

In ID the computations greatly simplify, since cells are 
mere intervals characterized by a single scale, their length 
I or equivalently their mass m = p§i. In agreement 
with Eq. (|70[) , the convex hull ip(q) between the boundary 
points qo and q± is the straight line 

<f(q) = VL;Q + — (PL;1 - <PL;0), (78) 

qi - qo 

the Eulerian position x of the cluster (shock) is 

fL;l - <PL;0 



<7i - 9o 



and its derivative reads as 
dx 



t 



d^o,o qi - qo 



(79) 



(80) 



k^T 1 : W (k£)~(k£) 



-d-l 



(74) 



The integral over s in Eq. (|57[) is now the discrete sum 
over the two boundary points, go and q\, and taking the 
average over the two orientations of the interval with re- 
spect to qo = 0, that is, either qi = £ or qi = —£, we 
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obtain 

f°° t 

G K j,(Aq,t) = / dmn(m)-[d D (Aq + £) + S D {Aq-£) 
Jo " 

(81) 
which 



-2S D (Aq)] . 

This is the result that was already obtained in [2c 
yields in Fourier space 

W(k;m) = [I-cob(W)], 



and 



G K (k,t) 



2t 

dmn(m) — [1 



(82) 



cos(fcf)] . (83) 



This also means that we recover the high- A: decay (|74|) 
and the high-A; power-law (|T6|) is indeed realized for all 
values of n in the range — 3 < n < 1 that we consider in 
this paper, because for high k the integral (|83| is always 
governed by the region £ ~ \/k. Thus, we recover the 
results 1(71)) and (fTl?)) that were obtained in Sect. MI C 31 
from simple dimensional and scaling arguments. This is 
not surprising since in ID Lagrangian cells are character- 
ized by a single scale. 



5. 2D case 

In 2D the computation follows the same route, the 
Lagrangian cells V being now triangular. Let us con- 
sider a triangle of summits {Si, S2, S3}, inner angles 
{ai, a-2, 03}, and opposite-side lengths {si, S2, S3}. We 
can also choose the triangle to be positively oriented, 
that is, S1S2 x S1S3 = S2S3 sin(ai)e3, where {ei, 62,63} 
is a 3D right-handed coordinate system and the triangle 
is in the plane {e!,e 2 }. Then, from Eq. (j70")) we obain 
the equation </?(q) of the planar facet that goes through 
the three summits, and from Eq. flTlT) the position of the 
mass cluster is x = dip/dq. Taking the derivative with 
respect to VtySu from Eq.©, we obtain 

<9x t ^ 

7T, = 5« S ^ x e 3' ( 84 ) 

dip - Sl 2£ 2 

where £ 2 is now the triangle area, and the scalar product 
with the outer normal n ou t along each triangle side reads 
as 



(SiS 2 ) 
(S2S3) 
(S3S1) 



{SxS 2 ) 


<9x 


—tS\ 


out 




2£ 2 


(S2S3) 


<9x 


ts\ 


out 




= 2P' 


(S 3 Si) 


9x 


—tS\ 


out 




2£ 2 



COS OL2, 



COS OL3. 



(85) 
(86) 
(87) 



Then, the kernel W reads from Eq. ([6H|) as 
3si 



W(k; m 



-si e 



2fc 2 

k-[qs 2 +"(qs 3 



J\u [S3 ooBoae-^* 



-1S 2 )] 



S2 COS Q!3 e 



-iukqs. 



(88) 
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FIG. 2: The resulting kernel function function W K {k;m) for 
equilateral triangles (dashed line) and for the angle distribu- 
tion of a Delaunay triangulation of a Poisson point process 
(solid line). 



where we chose the origin of coordinates as qsj = 0. Tak- 
ing the average over the global direction of the triangle 
(which is equivalent to averaging over the direction of k), 
we obtain 



W(k; m) 



3si 



dit 



S3 COS C*2 Jo(uks3) 



2k 2 £ 4 

S2 cos 0.3 Jo(wfcs 2 ) — Si 



Xj (fc v / S2 2 U 2 + S3 2 (1— U) 2 + 2S2S3U(1 — u) COSQfl) 



(89) 



We can be here a bit more explicit in the ensemble av- 
erage calculations. Let us define the angle distribution 
-Ptri.(o;i, ol 2 ', rn)da\da2 for a give mass (i.e. area) with 
the relations (with a3 = 7r — a\ — 02), 



I = —S2S3 sinai, 



s 2 = 2^ 2 (cot a x + cot a 2 ) 



(90) 



(91) 



and those obtained by circular permutations. As a result 
the expression (|59")l can be rewritten 



W(k;rn) = 



du / Ptn.{ai,a2;m)daida 2 



k 2 £ 2 J 

cot a 2 Jo(uks3) + cot CK3 Jo(uks2) — (cot a>2 + cot 0:3) 



X Jn(fc\/ S2 2 U 2 + S3 2 (1 — 7i) 2 + 2s2S3u(l — u) COStti) 



.(92) 



Note that expanding the Bessel functions up to order k 2 
and using standard relations between the angles, sides 
and area of triangles, the property (|6"T|) can be easily 
recovered. 

The actual distribution P t ri. to use actually depends on 
the details of the Lagrangian-space tessellation, which in 



11 



turn depend on the index n of the initial conditions. For 
instance, if all triangles were equilateral, we would have 



W e q(k;m) 



8 



k 2 s 2 



du [Jo(uks) — Jo(ksy/ u 2 — u + 1)] 

(93) 

with s 2 = (4to)/(v / 3/3o)- Another interesting, and proba- 
bly more realistic, case is provided by the angle distribu- 
tion obtained in the Delaunay triangulation of a Poisson 
point process, for which we have, (44. l45j. 

-Ptri.( a i: a2)daida 2 = 

■ sin(ai) sin(«2) sin(ai + a^daida^ (94) 



3tt 



with 



< ai < 7T, < tt2 < 7T, < U\ + OL1 < 7T, 



(95) 



independently of the area I 2 of the considered triangles. 
The resulting angle integration cannot be performed ex- 
plicitly. We instead propose a fit to the resulting dimen- 
sionless kernel, 



W(k;m) = 0.749428 e" ' 206885 ^ 2 
0.250572 



v/0.00707579(fcf) 6 + 1.00806(fc£) 2 + 1 



(96) 



The resulting shapes of the kernel are shown on Fig. [2] 

These forms illustrate the high k asymptotic forms of 
the propagators. At fixed triangle area and shape, the 
integration over the Bessel functions in (|89|) leads to a 
high-fc power-law decay of the form W(k;m) <~ (fc£) -3 , 
in agreement with Eq. (|74l) . This gives the n-dependent 
decay G K (k,t) ~ t(kL)- n -' 3 of Eq.([76]) for n < 0, and 
the constant-slope decay G K {k,t) ~ t(kL)~ 3 of Eq. ([77|) 
for n > 0. 

However, these exponents could be modified if the dis- 
tribution of triangle shapes, unlike the two cases con- 
sidered above, gives significant weight to configurations 
where the three sides have very different lengths. For 
instance, in the extreme case where triangles are increas- 
ingly "squeezed", with two sides that remain of order 
Lit) and a third side e that scales with mass as m ~ eL, 
the kernel is dominated by the case where the sum- 
mit 5*1 is at one end of the small side and it behaves at 
high k (and small mass) as 

m ~ eL, Ws qU cczcd(fc; m) ~ fc" 3 e" 2 £ _1 . (97) 

This now gives over the full range —3 < n < 1 the high-fc 
decay 

k » L^)- 1 : G K , squcczcd (fc, t) ~ t (kL)-^ n+ ^l\ (98) 



which is shallower than the scalings (|76| - (l77l) . This shows 
how the results of Sect. lIII C 31 which rely on simple scal- 
ing arguments and dimensional analysis, can be violated 
in complex cases where Lagrangian cells of a given mass 



involve several scales with different characteristic expo- 
nents. 

We can expect the scalings (f9"5]) and (|76l - ([77|) to 
bracket the high-fc exponents that can be reached for ar- 
bitrary triangle distributions, depending on scaling or on 
the fraction of "squeezed" triangles as a function of mass. 
This holds for mass functions that satisfy the low-mass 
power laws (|75p , but for more general cases we could ex- 
tend this analysis to arbitrary low-mass behaviors of the 
mass functions. In any case, these results show that the 
Lagrangian propagator is a sensitive probe of the struc- 
ture of the matter distribution. Moreover, in dimensions 
2 and higher, it involves both the cluster mass function 
and key properties of the shape of Lagrangian-space tes- 
sellations, that is, of the shape of the initial regions that 
eventually end up in small-mass clusters. 

Another significant difference with the Eulerian prop- 
agators discussed in Sect. IIIIBl which vanished for — 3 < 
n < — 1 because of infrared divergences of the initial ve- 
locity field, is that the Lagrangian propagators are finite 
and well-defined over the full range — 3 < n < 1. 



IV. NUMERICAL SIMULATIONS 

We now describe the results that we obtain from a nu- 
merical study, in both the ID and 2D cases. This allows 
us to check the exact results presented in Sect. IIIII and 
our predictions (for d = 2), which we think are illus- 
trative of what is happening at higher dimension. The 
details of our numerical simulations, and of the new effi- 
cient algorithms that we use, are given in |2l| . There, the 
statistical properties of the matter distributions built by 
the geometrical adhesion model defined in Sect. Ill Al for 
the same Gaussian scale- invariant initial conditions given 
in Sect. Ill Bl are described in details. In particular, it is 
shown that the density probability distributions and the 
mass functions exhibit qualitatively the same properties 
as those observed for 3D gravitational clustering (see also 

EE El)- 

Evolving from the initial Gaussian, that remains rele- 
vant on large scales and at early times, the density prob- 
ability distribution gradually broadens and builds an in- 
termediate power-law regime, in-between rare voids and 
rare overdensities. On the other hand, the shock mass 
function shows a low-mass power tail and a high-mass 
exponential-like falloff, and obeys up to a good accuracy 
the Press-Schechter-like scaling in terms of the reduced 
variable v = 5l/<j{M) [47j. However, while in ID the 
scaling mass function f(y) agrees quite well with the orig- 
inal Press-Schechter ansatz (which actually happens to 
be exact in the case n = — 2 [42]), in 2D it is significantly 
different and it shows a v 2 low-mass tail instead of the 
Press-Schechter prediction oc v. 

These studies have shown that the adhesion model 
shares many properties with the gravitational dynamics, 
which motivates a further investigation with respect to 
the response functions and correlators, which are criti- 
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cal ingredients in analytical approaches to gravitational 
clustering. 

The advantages for using the adhesion model are dis- 
cussed in [2~0 ] where one can further find the descrip- 
tion of the algorithms we use. Let us remind that our 
method allows us to get the mass distribution without 
introducing further discretization, either in time or in 
space. We can also cover a greater range of masses and 
scales as compared with usual gravitational dynamics (es- 
pecially as we consider lower-dimensional systems, cither 
ID or 2D), which gives us a better control of asymp- 
totic regimes. For numerical purposes, the power-law 
initial conditions (|20|) also improve the statistics since by 
using the self-similarity (|2Tj) we can rescale coordinates 
and propagators to obtain time-independent functions, 
so that we can take the mean over all output times when 
we compute statistical averages. 

More precisely, it is convenient to introduce the dimen- 
sionless scaling variables 

Q = mr * = Wr K = mK 

Then, equal-time statistical quantities (such as correla- 
tions or probability distributions) written in terms of 
these variables no longer depend on time and the scale 
X = 1 is the characteristic scale of the system, associated 
with the transition from the linear to nonlinear regime. 
In terms of the propagators, this gives for instance 

Gs(x,i) = -jj-r-j Gs{X), G s (k,t) = tG s (K). (100) 

In particular in the linear regime we have, 

Gf(X) = J D (X), Gf\K) = l. (101) 

The Lagrangian propagator G K obeys the same relations. 

As in [21j . we consider the cases n = 
0.5,0,-0.5,-1,-1.5,-2, and -2.5 to cover the 
range — 3 < n < 1 that we study in this article, where 
the self-similarity ([2~Tj) holds. This is also the range of 
interest for cosmological purposes. In particular, the 
definition of the slope n in Eq. (|2"0"|) is such that, whatever 
the dimension d that we consider, it corresponds to the 
usual slope n that appears in papers on 3D gravitational 
clustering in cosmology (thus, with this choice the 
scaling law (|2"Tj) does not depend on d). 

We focus on the Fourier-space propagators in the main 
text because they are the quantities that appear in the 
resummation schemes used in cosmology, However, we 
briefly describe in App. [B] our results in real space, for 
the response function G$ (X) . This is of interest by itself 
because for the Burgers dynamics this is also the one- 
point velocity probability distribution, as seen in Eq. (|4T|) . 
Moreover, in the nonlinear regime real-space methods 
may prove as useful as the Fourier-space methods that 
have been developed so far. 



A. Eulerian propagators 

1. Numerical estimates of propagators 

To compute the Fourier-space response function we use 
the definition (|2"5)l . as the mean functional derivative of 
the fields in Fourier space. Thus, taking as a reference 
an initial condition ^>o( x ), i.e. V>o(k) in Fourier space, 
we obtain the nonlinear field tp{x.,i), whence ip(k,t), at 
time t. Then, we perturb the initial condition at a given 
wavenumber ko . Since all fields are real, this means that 
we consider the change ipo(x.) — > ipo(x) + Ay»o(x) with 
Aipo (x) = e 2 cos(ko ■ x) , with a small amplitude e. Then, 
we compute the nonlinear field ip + Aip at time t gener- 
ated by this new initial condition, we take its Fourier 
transform, and we obtain the functional derivative as 
Dip/Vipo — Aip(ko)/e. Next, to estimate the response 
function G^(ko,t) as in Eq. ([2"5|) . with the statistical av- 
eraging (..), we repeat the operation above for many ref- 
erence initial conditions ipo, with the Gaussian statistics 
of Sect. Ill Bl and the mean over all these simulations gives 
our final estimate of G^(ko, t). In 2D we perturb the ini- 
tial condition along the two axis, ko = fcoei and next 
ko = ^0^2, which provides two measures of the func- 
tional derivative. Finally, using the self-similarity (|2"Tj) 
we rescale the results obtained at different output times 
and we make a final averaging to increase the statistical 
accuracy, to obtain the dimensionless propagators as in 
Eq. ([100j) . A similar procedure also provides the propa- 
gator Gs- Note finally that for the ID case and for the 
reduced variable K, we have 

G™(K) = G™(K), (102) 

but this identity is not true for higher dimensions, as 
shown in Sect. IIII Bl 

Instead of computing the response functions G(K) 
through functional derivatives as in Eq. (f2"51) , it is possible 
to derive these propagators through cross-correlations as 
in Eq. (|3T)]) . We also considered this alternative method 
to check our numerical algorithms. This also provides a 
confirmation of the general identities (|2"3"|) and (|27l) into 
the shell-crossing regime. We briefly show a comparison 
of these two procedures in App. [O for the ID case (we 
obtain similar but somewhat more noisy results in 2D). 

2. One- dimensional dynamics 

We show in Fig. [3] the Fourier-space Eulerian response 
function G^(K) = Gs{K), in terms of the reduced 
wavenumber K of Eq. ([Ml) for the ID case. From Eqs. fiTj) 
and (|48p. we have in real space 

G S (X) = G^(X) = P(U), with U = X, (103) 

where P{U) is the one-point probability distribution of 
the reduced velocity U, so that Gs(K) and G^(K) are 
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FIG. 3: The ID Fourier-space propagator Gg(K) = G^(K) 
from numerical simulations. We show the cases n = 0.5, 0, 
and —0.5, as a function of the reduced wavenumber K. For 
n — we also plot the exact analytical result (dot-dashed 
line) from [13]. 



also the Fourier transforms of the velocity distribution. 
Since for — 3 < n < — 1 the response functions are not 
well defined because of the divergent "sweeping effect" , 
we only show our results for the cases n = 0.5,0, and 
—0.5. For the case n = we also plot the exact re- 
sult, given by Eq.(65) and Figs. 1 and 2 in (23|. We can 
check that our numerical result agrees quite well with 
this analytical result and its strong exponential-like de- 
cay ~ e~ K 1 . However, the numerical errorbars are too 
large to distinguish the oscillatory behavior (i.e., changes 
of sign) in the far tail, where G < 10~ 3 . Our numerical 
results show that this response function obeys a simi- 
lar exponential-like decay at high K for n = 0.5 and 
n = —0.5. A priori the falloff can depend on n, as 
G ~ e~ K with some exponent a(n). The range of scales 
available in Fig.[3]is too small to obtain a precise measure 
of a(n), but it suggests that a does not vary too much 
over —0.5 < n < 0.5. 

For K — > we recover G^ — > 1, in agreement with 
Eq. (|4"o]) and with the normalization to unity of P{U) in 
Eq.CD!]). 



3. Two-dimensional dynamics 

Let us now consider the 2D case. We show in Fig. [4] the 
2D Fourier-space propagators Gs(K) and G^{K). As in 
the ID case of Fig. [31 we can clearly see the steep expo- 
nential cutoff, but with a slightly stronger dependence on 
the exponent n of the initial conditions. We also recover 
the exact low- if limit G^(0) = 1. 

As explained in Sect. IIIIB1 in dimension greater than 
one the density and velocity potential propagators are 
no longer identical, because of the nonlinear relationship 
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FIG. 4: The 2D Fourier-space propagators Gs(K) (solid line) 
and G^,(K) (dashed line). 

(f5U|) . but they are all expected to be governed by the 
"sweeping effect" . This is confirmed by our numerical re- 
sults, since we clearly see in Fig.|4]that Gs(K) shows the 
same exponential-like cutoff at high K. In fact, in 2D the 
density propagator is very close to its velocity-potential 
counterpart, although there are hints of a slightly steeper 
falloff for the density response function. 



B. Lagrangian propagators 

As for the Eulerian propagators studied in Sect. II V ~K[ 
we measure the response function G K {K) from its defi- 
nition p8p. by estimating the functional derivative from 
the difference between two very close initial conditions. 
In contrast to the Eulerian propagators, the Lagrangian 
propagators are well defined over the full range —3 < n < 
1, hence we plot the cases n = 0.5, 0, —0.5, —1, —1.5, —2, 
and —2.5. 



1. One- dimensional dynamics 

We show in Fig. [5] the Fourier-space Lagrangian prop- 
agator G K (K). For the two cases n — and n = — 2 we 
obtain a very good agreement between our numerical re- 
sults and the exact analytical results derived in [23[ (the 
numerical and analytical curves cannot even be distin- 
guished in this figure). We recover the exact low- if limit 
G(0) = 1, which corresponds to Eq. (|rj5]) . As explained 
in Sect. IIII C 21 this is related to the relative motions of 
well-separated regions by long-wavelength modes of the 
velocity field, and it applies to the full range —3 < n < 1. 

Figure [5] clearly shows a power-law tail at high K, 
contrary to the exponential-like cutoff seen in Fig. |31 
with a slope that depends on n. To clearly see this 
power-law tail we show in Fig. [B]the logarithmic deriva- 
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FIG. 5: The Fourier-space response function G K (K) (solid 
line) from ID numerical simulations. For n — and n — — 2 
we also plot the exact analytical result (dot-dashed line) from 
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FIG. 6: The logarithmic derivative, — din G re /d In K, of the 
Lagrangian propagator, from ID numerical simulations (solid 
line), and its theoretical asymptotic limit (n + 3)/2 (dashed 
line). 



FIG. 7: The 2D Fourier-space propagator Gk(K) (solid line) 
from numerical simulations. 



... 




FIG. 8: The logarithmic derivative, — dlnd/dln K, of the 
Lagrangian propagator, from 2D numerical simulations (solid 
line). We also show the two theoretical predictions of Eq. (|105p 
(upper dashed line) and Eg . (I106fl (lower dot-dashed line). 



2. Two-dimensional dynamics 



tive — d In G K /d In K, as well as its asymptotic theoretical 
prediction from Eg. (|83|) . which agrees with Eq. ([76| and 
reads as 



d=l, -3 < n< 1, K -> 



dlnG K 
din A" 



(104) 

We can check that our numerical results are consistent 
with this analytical result and we clearly see the conver- 
gence to the high- A" slope (| 104[) . 



We show in Fig. [7] the Fourier-space Lagrangian re- 
sponse function G K (K). As in the ID case plotted in 
Fig. [5J we recover the low- AT limit G K (0) = 1. We again 
obtain a clear power-law tail at high AT, but the conver- 
gence to this asymptotic regime is somewhat slower than 
in ID, especially for low values of n. 

To estimate the high- AT exponent we plot in Fig. [5] the 
logarithmic derivative, — dlnG K /dlnAT. As explained in 
Sect. IIII C 51 depending on the shape of the Lagrangian- 
space tessellation, we have two different predictions for 
the extreme cases of small-mass triangles characterized 
by a single scale I ~ to 1 / 2 , that is, triangles of approxi- 
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mately "equilateral" shape, 



"equil." , K — > oo 



dhxG K 
din if 



-min(n + 3,3), (105) 



which follows from Eas. ([75|) . (|77)) . and for "squeezed" 
triangles, where two sides remain on the order of the 
nonlinear scale L(t) whereas the third side decreases as 



"squeezed" , if 



din G K 
din if 



3(n + 3) 



(106) 



from Eq. tjM)) . We clearly see in Fig. [5] the convergence 
towards a constant logarithmic slope, although for low n 
the convergence is quite slow and is not complete yet on 
this range. As expected the high- if exponent is between 
the two limiting values (|105l) and (|106l) . Unfortunately, it 
also seems that it is different from these two values, which 
means that the triangulation is neither dominated by 
"equilateral" (i.e. "single-scale") triangles nor by "max- 
imally squeezed" triangles. The Lagrangian tessellation 
is probably more complex than these two simple cases 
and Fig. [5] suggests that there is a broad distribution of 
shapes for low-mass triangles. 

This agrees with a qualitative inspection of Fig. 6 in 
[2(1. where we showed two examples of Lagrangian tri- 
angulations obtained for n = and n — —2. There, one 
can clearly see that the summits of the Lagrangian cells 
are not distributed at random but show a strong cluster- 
ing. One can see distant groups of Lagrangian summits, 
separated by a length of order L(t), with small triangles 
within each group (they would correspond to low-mass 
"equilateral" shapes) and very thin triangles that join 
two summits in one group to a third summit in a second 
group (they would correspond to low-mass "squeezed" 
shapes) . This suggests that the Lagrangian-space tessel- 
lation contains both "equilateral" and "squeezed" con- 
figurations. However, it is not clear whether there is a 
bimodal or a continuous distribution of triangle shapes. 

For completeness, we should point out that numerical 
results suggest that the number of mass clusters per unit 
Lagrangian or Eulerian volume is infinite if— 3<n<— 1 
and finite if — 1 < n < 1, as seen from the low- mass tail 
(l75j) of the mass function [13, [HI- Moreover, it appears 
that clusters are dense in Eulerian space for — 3 < n < — 1 
while they are isolated for — 1 < n < 1. In Lagrangian 
space the geometry is somewhat more complex. For 
— 1 < n < 1 there is still a finite number of cells and 
summits per unit area, but for — 3 < n < — 1 the infi- 
nite number of triangles and summits does not cover the 
whole plane. This is obvious from the fact that there 
are some large triangles associated with finite-mass clus- 
ters. Therefore, the infinitely many summits "gather" in 
some regions of Lagrangian space, in-between the large 
triangles associated with massive halos[55f. Nevertheless, 
these rather different properties of the Eulerian and La- 
grangian tessellations with n [2(| do not seem to have 
a strong impact on the low-mass tail nor on the high- 
if asymptote of the Lagrangian response function (e.g., 
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FIG. 9: The statistical average (L\) of the longest triangle 
side L\. We plot this quantity for a series of triangle mass 
bins and we divide {Li) by the factor yM, where M is the 
mean mass of each bin. 




FIG. 10: The statistical average (L3} of the shortest triangle 
side L3. We plot this quantity for a series of triangle mass 
bins and we divide (L3} by the factor \J~M, where M is the 
mean mass of each bin. 



there is no clear sign in Fig. [8] that the exponent switches 
from (I105P to (|106l) as n goes to either side of —1). 



3. 2D Triangle shape distribution 

As discussed in the previous section an important in- 
gredient for the derivation of the Lagrangian propagator 
is the triangle shape distribution function. In particular 
the derivation of Eq. (1105)) assumes that the triangles that 
build the Lag rangian-space tessellation of the matter dis- 
tribution [2QJ are scale-invariant in the small-mass limit. 
In other words, it neglects any dependence of their shape 
probability distribution on scale since it assumes that 
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typical lengths scale as the square-root £ of the triangle 
area (in the limit of small triangles) . This assumption is 
violated if the shape distribution does not converge to a 
finite limit, for instance if smaller triangles are increas- 
ingly "squeezed" as in Eq. (|106p . 

To investigate this point, we consider in this section 
the dependence on scale of the triangle geometry. Thus, 
for each realization and at each output time, we label 
the three sides of each triangle of the Lagrangian-space 
tessellation according to 

ti>ii>ia, (107) 

that is, l\ is the longest side and £3 the shortest side. 
Then, we plot in Figs . l9l and ITUl the averages of the ratios 
Li/VM = h/VA and L 3 /VM = £ 3 /VA, where A = £ 2 
is the triangle area. To do so, we first bin the dimen- 
sionless triangle area M (which is also the dimensionless 
mass) over a finite number of bins. Then, for each re- 
alization we compute the dimensionless lengths L\ and 
L3 and mass M of all triangles found in the Lagrangian- 
space tessellation at a given time, we count all triangles 
that fall within a given mass bin and we also store their 
lengths L\ and L 3 . Repeating this operation over sev- 
eral output times (taking advantage of the self-similarity 
(|2ip ) and over many initial conditions, characterized by 
the Gaussian statistics described in Sect. Ill 51 we obtain 
the means (L\) and (L 3 ) within each mass bin, whence 
the ratios (Li)/y/M and (L 3 )/y/M plotted in Figs.[3]and 

EES 

If the triangle distribution is scale-invariant in the 
small-mass limit (i.e. for small area) these ratios must go 
to finite and nonzero constants at low M. Figures |H] and 
[10] suggest that this is not the case, and that (Li)/y/M 
keeps increasing while (L 3 ) /vM keeps decreasing at low 
mass. This means that triangles become increasingly 
"squeezed", that is, £\ ~ £ 2 and £3 <^ £1. This behavior 
is consistent with a visual inspection of Fig. 6 in [201 . 
It also agrees with the analytical result obtained in |48j 
when the velocity potential tJjq (q) is a random Poisson 
point potential. However, this corresponds to very differ- 
ent initial conditions than the ones studied in this paper, 
since then tpo(q) is drawn from a finite probability dis- 
tribution at each point and there are no correlations be- 
tween different points. In contrast, in our case there are 
significant large-scale correlations, which are character- 
ized by the exponent n. As discussed in Sect. lIVB 21 this 
leads to strong correlations in the Lagrangian-space tes- 
sellation itself, with a combination of "equilateral" and 
"squeezed" triangles that join summits that are in the 
same or two well-separated groups. 

We can see in Figs . [91 and [TU1 that the scale-dependence 
of the mean triangle shape is however quite weak and 
only logarithmic (or possibly a power law with a small 
exponent), since the ratio {Lx)/yM typically grows by a 
factor 5 and (£3) /vl decreases by a factor 2 when M 
decreases by 4 orders of magnitude. This would validate 
the scaling prediction (|105[) . up to logarithmic predic- 
tions. However, the deviations seen in Fig. [5] suggest 



that the mean values (Li) and (L 3 ) are not sufficient to 
derive the precise behavior of G K , which is quite sensitive 
to the relative fractions and "squeezed" and "equilateral" 
configurations, and of the intermediate shapes, since dif- 
ferent exponents for L3 as a function of M give rise to dif- 
ferent exponents for the high-X tail, in-between the two 
limiting cases (I105P and (|106[) . Nevertheless, the slightly 
steeper dependence on mass of the ratios (Li)/yM and 
(L 3 ) /VM observed for larger n agrees with Fig. [8j where 
the asymptotic slope of the propagator appears to move 
farther from the "equilateral" prediction for larger n. 

Here we should note that these results differ from the 
properties associated with the Delaunay triangulation of 
a Poisson point process, which are recalled in App. [A] 
This is not surprising, since the Lagrangian-space tesse- 
lation is not a Delaunay triangulation and the supporting 
points are non-Poisson distributed. Thus, this discrep- 
ancy is another signature of the non-negligible correla- 
tions that are the results of the power-law initial condi- 
tions (I2TJ1) . which provide important large-scale correla- 
tions. 



V. DISCUSSION AND CONCLUSION 

In this paper we have explored the concept of propaga- 
tors for the geometrical adhesion model, paying special 
attention to the Lagrangian propagators. As recalled in 
the introduction propagators are defined in general as re- 
sponse functions of the final density field, or velocity field, 
to an infinitesimal change of the initial conditions. For 
Gaussian initial conditions these propagators can also be 
expressed in terms of unequal-time correlation functions. 
Following [13] we show in the first Section that this iden- 
tity is quite general, more general than what was initially 
thought in |6|, and can be derived beyond perturbation 
theory calculations. 

The main result of this paper is the relation between 
the Lagrangian propagators and the halo mass function 
that we uncovered in the context of the GAM. The key 
point with which we established this connection is that, 
as soon as all the particles are gathered in halos, the dis- 
placement field in Lagrangian coordinates describes the 
formation of a partition where each cell correponds to 
a single halo. Then, the convex hull construction that 
underlies the GAM explicitly shows that an infinitesi- 
mal variation of the initial conditions does not induce 
a change in this Lagrangian partition, or in the mass 
of each halo, but only modifies the halo final positions. 
This can be understood by noticing that propagators are 
defined as "linear" response functions, in the sense that 
they describe the sensitivity of the system to perturba- 
tions up to linear order (but include the full nonlinear 
background dynamics). Although we did not do it ex- 
plicitly it is clear that they can also be generalized to 
higher orders by expanding the response over powers of 
the perturbation. It is also clear that this result is not 
specific to this model. It should be valid as soon as the 
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dynamics results in the partition of the Lagrangian space 
into cells associated with distinct halos. 

With the key property described above, the functional 
relation (joTj) between the halo mass function and the 
propagator takes the form, 



G K (k) 



dmn(m) m W(k; m) 



(108) 



This provides an explicit link between the low-mass tail 
of the halo mass function and the high-fc tail of the 
propagator. Moreover, the dimensionless kernel W(k;m) 
contains some further information on the low-mass La- 
grangian cells, more precisely it describes the scalings 
with mass of their typical lengths (in dimension greater 
than one there may be more than one relevant scale, as 
cells may become more or less "squeezed" or "flattened" ) . 
This relationship provides a well-defined range of possible 
exponents for the asymptotic behavior of the Lagrangian 
propagator, if we know the low-mass tail of the halo mass 
function, and a definite prediction if we also know the 
shape of low-mass cells. 

In a more general context, such as the 3D gravita- 
tional dynamics, we cannot derive an explicit relation- 
ship of this form. However, it is natural to expect again a 
strong link between the Lagrangian propagators and the 
halo mass function. This must be contrasted with the 
Eulerian propagators, where in both dynamics the main 
process at play for these initial conditions is a "sweeping 
effect" due to long-wavelenght modes of the velocity field. 
Thus, our results strongly suggest that Lagrangian prop- 
agators are generally much better probes of the density 
field, and in particular of the halo mass function. More- 
over, a new result that we have uncovered in this paper 
and that appears in dimensions greater than one is the 
sensitivity to the shape of the underlying Lagrangian- 
space tessellation, through the kernel W(k;m), which 
should also remain valid for more general cases. 

It would be interesting to use the deeper understand- 
ing brought by such studies to improve or build quantita- 
tive tools. With respect to the Lagrangian propagators 
two avenues naturally appear. One could first use La- 
grangian response functions to probe some properties of 
the system. However, in a cosmological context where 
one cannot perform actual experiments (except for nu- 
merical simulations) this may not be very practical. A 
second route would be to use our understanding to build 
more efficient analytical schemes, such as new resumma- 
tion methods of Lagrangian perturbation theory. 



Appendix A: Angle distribution in Delaunay 
triangulation from a Poisson point process 

We consider a set of points with a Poisson distribution 
in a d = 2 space and the Delaunay triangulation associ- 
ated with those points. It is then possible to derive the 
joint probability distribution function of the three angles, 



ai , a? and ai3 

by mm, 



7r — a% — a2, of the triangles. It is given 



P t ri.(ai,a2)daida 2 = 
g 

— sin(ai) sin(o!2) sin(ai + a2)daida2 (Al) 

07T 



with 



< ttl < 7T, < Q!2 < 7T, < ttl + «2 < 7T. 



(A2) 



This distribution entirely characterizes the shape distri- 
bution of the triangles. In particular one can derive the 
average length of each side for a given area. Each side 
length is given by 



A 



2 sin 
sinai sin ai 



(A3) 



where A is the surface area of the triangle and where the 2 
other side lengths are obtained by circular permutations 
of the indices. The resulting average value of the longest 
and shortest sides are then respectively, 



(Li) 



A 



and 



(L 3 



A 



2.39298 



1.12187. 



(A4) 



(A5) 



These results are to be compared to those obtained in 
our construction. 



Appendix B: Real-space Eulerian propagators 

As seen in Eq. pTt the real-space response function 
G^,(x, t; qo) of the velocity potential is given by the one- 
point velocity probability distribution function. In terms 
of the reduced variables (|99|) . using the statistical homo- 
geneity of the system, this reads as 



GV(X) = P(U), with U = X, 



(Bl) 



which holds in any dimension. We show in this appendix 
our results in ID and 2D for this propagator, or veloc- 
ity distribution, which is also the inverse Fourier trans- 
form of the propagators plotted in Figs. [3] and HI To 
measure this real-space quantity we do not use the func- 
tional derivative (|3"TT) . Instead, using Eq. (|Bip we simply 
measure the velocity distribution. 

We can note that the identity (|B1[) implies that the 
real-space propagator G^ is always positive. This is not 
the case for its Fourier-space counterpart G.0. For in- 
stance, it was proved in 23j that in ID for n = the far 
tail of G^(K) shows fast oscillations that are exponen- 
tially damped (a first oscillation can be distinguished on 
the dot-dashed curve in Fig. [3]) . 

As in Sect. II V Al we only consider the initial conditions 
n = 0.5,0, and —0.5, since as explained in Sect. MI Bl 
these Eulerian propagators and one-point velocity distri- 
butions do not exist for n < — 1. 
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FIG. 11: The real-space response function G^,(X) = Gg(X) 
from ID numerical simulations (solid line). We show the cases 
n = 0.5, 0, and —0.5, as a function of the reduced distance X. 
For n = we also plot the exact analytical result (dot-dashed 
line) HI]. 




X 

FIG. 12: The ratio - ln{G i ,(X)]/X (n+3 '> from ID numerical 
simulations (solid line). For n = we also plot the exact 
analytical result (dot-dashed line) and its asymptotic limit 
(dashed line) [H]. 



1. ID case 

We show in Fig. [TT] the real-space response function 
G^(X) = G S (X). As in Sect. IIV A 2l we can check that 
for n = our numerical result agrees very well with the 
exact result (which cannot even be distinguished in this 
figure), given in (2^,|H,|3^. The response functions ob- 
tained for n = 0.5 and —0.5 show the same behavior, 
with a sharp large-distance falloff, but contrary to its 
Fourier transform plotted in Fig. [3] we can clearly see the 
dependence on n (especially for n = —0.5). 

To clearly see this large-distance falloff, we plot in 
Fig. [Hthe ratio - \n[G 4 ,{X)]/X < ~ n + 3 \ which must go 



FIG. 13: The real-space response function G,/,(X) = P(U), 
with X = U. We show the probability distribution of each 
velocity component, P(Ui) (solid line), and the probability 
distribution of the total amplitude, P(|U|) (dashed line). 



to a constant at large X from the analytical result (|44p . 
As in Fig. 111! we can check that for n = our numeri- 
cal result follows the exact result from [2^, [3^, HH . The 
curves obtained for n = 0.5 and —0.5 display a similar 
behavior, but as shown by the explicit case n = the 
curves have not converged to their asymptotic limit yet. 
Therefore, we cannot precisely measure from the simula- 
tions the exponent of the large-distance falloff, but they 
are consistent with the analytical result (|4"4"]l . The values 
reached in Fig. [T2l correspond to G^{X) and P(U) below 
10~ 10 for the case n — 0, which means that the conver- 
gence of the ratio — \n[G^(X)]/X^ 1+3 ^ to its asymptotic 
limit is rather slow. 



2. 2D case 

We show in Fig. [13] the real-space response function 
G^ (X) , or more precisely the velocity probability distri- 
bution P(U), using the identity (IB1|) . This is the in- 
verse Fourier transform of the response function G^(K) 
shown in Fig. [4] As in the ID case, we obtain a sharp 
exponential-like cutoff at large distance (i.e. at large 
velocity in terms of P(U)), which is consistent with 
Eq.flUD. 



Appendix C: Cross-correlations 

We show in Figs. [T4l and [T5l the relative difference be- 
tween the response functions G(K), computed through 
functional derivatives as in Eq. (|28p . and the same quan- 
tities computed through cross-correlations as in Eq. (f3"0"f , 
which we denote as S(K) to distinguish between both 
procedures. We show our results for the ID Eulerian 
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K 



propagator in Fig. Qj] and for the ID Lagrangian 
propagator G K in Fig. \T5\ We plot the upper limit that 
we obtain for the relative difference \G— S\/\G\, since our 
errobars are everywhere consistent with the exact equal- 
ity G = S. (Of course this upper limit weakens at high 
K where the propagators are very small and numerical 
measures are more difficult.) Wo obtain similar results 
in 2D but the tests of this identity are somewhat weaker 
because of the lower-quality statistics. 



FIG. 14: The upper limit for the relative difference \G$ — 
S^,\/\G^,\, from ID numerical simulations. We show the cases 
n = 0.5, 0, and —0.5, as in Fig. [3] 




FIG. 15: The upper limit for the relative difference \G K 
S K \/\G K \, from ID numerical simulations. 



These figures show the validity of our numerical algo- 
rithms and provide an estimate of the accuracy of our 
measures of the propagators. They also confirm the va- 
lidity of the identities (|2"5)) and (|27l) in the nonpertur- 
bative highly nonlinear regime, in agreement with the 
derivation presented in the introduction. 
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Here "linear" refers to the fact that we only consider the 
response up to linear order over i.e. for infinitesimal 
external perturbations. However, Gs( includes the fully 
nonlinear dynamics of the physical field S itself, 
at least in the case of single pressureless fluid. This is not 
the case for multiple fluids, see [32T ]. 

Here we used the standard definition of the Legendre- 
Fenchel conjugate f*{s) of a function /(x), /*(s) = 
£ s [/(x)]=sup x [s-x-/(x)]. 

This amounts to linearize the equation of motion (|3} and 
the continuity equation (which holds before the formation 
of shocks). 

This only holds in the range —3 < n < 1. 
As noticed in Sect. MI A 31 the divergence k is not sensi- 
tive to uniform translations of the system, hence it may 
seem puzzling that the normalization (|69[) involves the 
Galilean transformation (|64p . However, this can be un- 
derstood as follows. The large-scale limit (|69[) an d the 
linear behavior (|34jl are direct consequences of the re- 
quirement that two far-away regions xi and X2 , with ini- 
tial velocities uo ; i and uo ; 2, see their relative distance 
evolve as X2 - xi = i (uo : 2 — u o,i). This sets the ampli- 
fication factor t in Eq. (|34[) . This configuration considers 
that regions 1 and 2 have constant velocities on some 
neighborhood and neglects higher order effects, so that 
they move independently according to the Galilean law 

dSl. 

These Eulerian- and Lagrangian-space behaviors are also 
found in the ID case, where they are rigorously proved 
for n — —2 and n = 0. 



